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Abstract 

A method is developed for estimating the emission rates of contaminants into the atmosphere from multiple point 
sources using measurements of particulate material deposited at ground level. The approach is based on a Gaussian 
plume type solution for the advection-diffusion equation with ground-level deposition and given emission sources. 
This solution to the forward problem is incorporated into an inverse algorithm for estimating the emission rates by 
means of a linear least squares approach. The results are validated using measured deposition and meteorological data 
from a large lead-zinc smelting operation in Trail, British Columbia. The algorithm is demonstrated to be robust and 
capable of generating reasonably accurate estimates of total contaminant emissions over the relatively short distances 
of interest in this study. 
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1. Introduction 

Urban air quality is an issue of major concern ow- 
ing to recent upward trends in population growth 
and urbanisation and industrialisation around the 
world. Consequently, there is an increasing need to 
understand the detailed dynamics governing emis- 
sion and transport of particulate matter in the at- 
mosphere. Turner [IJ mentions a multitude of pos- 
sible sources of air-borne particles, including those 
of anthropogenic origin such as industrial complexes 
and automobiles, as well as natural sources such as 
dust storms and volcanic eruptions. Recently, there 
has been a surge of interest in related problems 
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for disaster planning and national security that in- 
volve transport of radionucleotides and biological or 
chemical agents 

The physics of particulate transport in the atmo- 
sphere are complex, in many cases involving mul- 
tiple spatial scales (ranging from the particle scale 
to near-source and long-range effects), multi-physics 
(coupling mass transport, turbulence, chemistry and 
wet/dry deposition), and complex geometry (e.g., 
involving flow over topography or man-made struc- 
tures). Models for these and other aspects of atmo- 
spheric dispersion have a long history dating back 
to the pioneering studies of turbulent diffusion by 
Richardson [5] and Taylor [3] . The bulk of previous 
work has tackled the forward problem^ by which we 
refer to the process of determining downwind con- 
taminant concentrations given source emission rates 
and meteorological conditions. These forward mod- 
els are usually based on a solution of the advection- 
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diffusion equation that is obtained through either 
analytical or numerical means (or a combination of 
both). The most prevalent approach used in prac- 
tice, and which is implemented in many industry- 
standard software packages, employs an approxi- 
mate analytical solution for point source emissions 
known as the "Gaussian plume solution." The one- 
dimensional plume solution for a single point source 
was originally derived by Sutton |^ and has since 
been extended to higher dimensions, as well as be- 
ing applied to a variety of other more general situa- 
tions involving ground- level deposition 6J , multiple 
sources [7] , height-dependent wind speed and diffu- 
sion coefficients |8I9) . and line and area sources, to 
name just a few. We remark that while a major- 
ity of the applications considered to date have in- 
volved transport in the atmosphere, the Gaussian 
plume models just mentioned may also be used to 
solve other advection-diffusion problems in such di- 
verse areas as population growth [lOj or water flow 
in rivers pTj and the subsurface [1^ . 

Another related stream of research has focused on 
solving the corresponding inverse problem, whereby 
measurements of particulate concentrations or 
ground-level depositions are given and the aim is 
to determine information about the location or ef- 
flux rate of contaminant sources. Inverse methods 
based on Gaussian plume type solutions have been 
developed by a number of authors in this context 
including Jeong et al. Jl3j and Hogan et al. [14], 
while MacKay et al. [I5j developed an alternate 
solution approach using complex variable theory. 
Other researchers have applied a more direct com- 
putational approach by solving the nonlinear gov- 
erning equations using methods based on Kalman 
filtering [16j . Lagrangian particles [T7j . Bayesian 
techniques [18|19j . or by integrating the equations 
backward in time |17|20j . Several related meth- 
ods have been developed specifically for handling 
the added nonlinearities arising from atmospheric 
chemistry, typically using Newton type iterative 
methods [21j, and sometimes combined with sta- 
tistical techniques [22] • These direct numerical 
approaches can be very computationally intensive, 
especially for 3D problems, and so will typically 
require use of parallel computing resources. Re- 
gardless of the numerical method used, the inverse 
problem is often characterized as ill-conditioned in 
the sense that small changes in parameters can lead 
to very large changes in emission estimates; these 
issues are discussed in much more detail by Enting 
[18], Beychok [23], and El Badia et al. [TT] . 



The subject of this paper is a lead-zinc smelt- 
ing operation located in Trail, British Columbia, 
Canada. We are concerned with the transport of sev- 
eral contaminant species from multiple point sources 
on the site, and our aim is to develop an inverse 
algorithm that will determine emissions based on 
the Gaussian plume solution of Ermak [5] . The al- 
gorithm should be capable of generating reliable es- 
timates of emissions rates given a relatively small 
number of deposition measurements. The novelty of 
this work stems from a combination of factors: 

- We make use of real (noisy) meteorological and 
deposition data, in contrast with some other stud- 
ies that use synthetic data [14115] . 

- Deposition measurements are relatively small in 
number and represent time-averaged accumula- 
tions. This should be compared with some other 
methods that obtain very high accuracy by using 
very large numbers of sample points [13] . Another 
example is the work of Hogan et al. [14J, who pro- 
posed an iterative method based on the Gaussian 
plume solution (with constant wind and no depo- 
sition) and which exploits the fact that concen- 
tration measurements at four locations uniquely 
determine the location and strength of a single 
point source. This approach can be very effective 
when the input data are known very accurately, 
but it degrades when the data are noisy. 

- The emission sources are at known locations, in 
comparison with some other studies that aim 
to determine both emission rates and locations 
|l(il2()l21j . 

- We incorporate additional linear constraints on 
emission rates that are derived from chemical pro- 
cesses within the smelting operation. 

- Deposition measurements are taken near ground 
level and at short distances from the source, which 
allows us to avoid errors inherent in long-range 
dispersion estimation and thereby minimize the 
ill-conditioning of the inverse problem. 

Taken together, these factors allow us to develop a 
robust algorithm that is capable of estimating emis- 
sion sources with a reasonable degree of accuracy. 
Other studies have been performed on emissions at 
the Trail site by Goodarzi et al.(see [24] and refer- 
ences therein), but they use a much simpler Gaus- 
sian plume solution with no deposition and constant 
wind velocity, as well as validating their results us- 
ing long-range deposition measurements and differ- 
ent experimental techniques. 

We begin in Section [2] by describing the prob- 
lem under study and developing a detailed list of 
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assumptions underlying the model. In Section [3] 
we provide details of the Ermak solution to the 
advection-difFusion equation, and also incorporate 
multiple sources and a time-varying wind velocity. 
Section |4] focuses on the inverse problem, deriv- 
ing the linear equality and inequality constraints 
and describing the linear least squares solution 
algorithm. A scries of numerical simulations are 
performed in Section [5l including a study of the 
sensitivity of the model to changes in parameters 
and noise in the data. Finally, we conclude with 
recommendations about the suitability of applying 
our model to actual environment reporting scenar- 
ios, and make suggestions on possible future work 
on extending our approach in order to improve ac- 
curacy and permit application to a wider range of 
atmospheric dispersion scenarios. 

2. Problem description and simplifying 
assumptions 

The motivation for this work was a study of 
emissions from a number of contaminant sources at 
a large lead-zinc smelter located in Trail, British 
Columbia, Canada and operated by Teck Cominco 
Limited. Our primary aim was to improve the ac- 
curacy of airborne emission estimates (especially 
for zinc) that the Company is required to report 
annually to Environment Canada's National Pollu- 
tant Release Inventory (NPRI) [25j . Stack emissions 
on the smelter site are directly measured for zinc 
and other contaminants, but this paper deals with 
low-level sources for which direct measurements are 
difficult to obtain. 

There are four sources on the Trail site that emit 
zinc (in the form of zinc sulphate, ZnS04) and these 
are indicated on the aerial photo in Fig. [1] by the 
symbol Ss, where s = 1, 2, 3, 4. To assist in estimat- 
ing the level of zinc emissions, the Company per- 
formed a series of ground- level measurements of zinc 
as well as a number of other contaminant species 
(strontium, sulphur, etc.). The measurements were 
taken over the two-year span 2001-2002 and consist 
of one-month accumulations of particulates within 
dustfall jars or "receptors," which are located at nine 
separate locations Rr, r = 1, 2, . . . , 9 (also indicated 
in Fig. [J). 

Meteorological data is available for the same 
monthly periods in terms of wind speed and direc- 
tion averaged over 10-minutc intervals. The smelter 
is located in the Columbia River valley which tends 



to funnel the winds on the site in a specific direction; 
since the river is located just below the aerial photo 
in Fig. [1] and runs roughly horizontally, the prevail- 
ing wind directions are roughly to the northwest or 
southeast. 

Based on the above description, we now make a 
number of assumptions which are critical in the de- 
velopment of our dispersion model: 

(i) Each emission source is considered to be a 
point source, and all emission rates are con- 
stant in time, at least for any one-month pe- 
riod over which depositions are measured. 

(ii) The wind velocity u{t) depends on time only 
and is uniform throughout the domain at any 
given instant. This is justifiable considering 
the relatively small dimensions of the site 
(1600 X 800 m) and the short time intervals 
of interest. 

(iii) Variations in topography are negligible so that 
the wind can be assumed horizontal. Although 
Trail is located in a river valley bordered by 
steep mountains, the domain of interest is far 
enough removed from the neighbouring moun- 
tain range that any topography or boundary 
effects can be ignored. 

(iv) A ten-minute averaging period (or time step) 
is used in all calculations, which is consistent 
with the assumptions employed in deriving the 
Gaussian plume solution |23I26| to limit errors 
in concentration. 

(v) Only dry deposition is considered since the 
dustfall measurements were all taken during 
months for which rainfall is relatively small. 

(vi) The effects of plume rise are incorporated by 
using an "effective height" for each stack. 

(vii) The atmospheric stability class is assumed 
to be class D ("neutral", according to the 
Pasquill-Gifford classification scheme) for all 
monthly periods. There are insufficient me- 
teorological data available to consider vary- 
ing the stability class with time, and so our 
choice of neutral class is a compromise that 
takes into account predominant atmospheric 
conditions during the months of interest. 

It is important to note that even though we focus 
our attention on an application to a specific indus- 
trial site, the mathematical model and associated 
numerical algorithms we develop are general and so 
can potentially be applied to a wide range of other 
atmospheric dispersion problems. 
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Fig. 1. An aerial photo of the Trail site, showing the approximate locations of each of the sources (red circles, labeled by Ss, 
= 1,2,3,4) and receptors (green triangles, labeled Rr, r = 1, 2, . . . , 9). The size of the area depicted here is approximately 
1600 X 800 m and the directions of both the prevailing winds and compass north are indicated in the lower right corner. 



3. Derivation of the forward model 

The release and transport of a single contami- 
nant in the atmosphere can be described by the 
advection-difFusion equation 

ct + M • Vc = V • (/CVc) + S, (1) 

where 

c{x,t) = contaminant concentration [kgm^"^], 

u — convective wind velocity [ms~^], 

JC{x) — diag(i^i,, Ky, K^), the matrix of tm-bu- 
lent eddy diffusivities [m^ s^^], 

S{x, t) — emission source term [kgm^"^ s^^]. 
We begin by considering a single elevated point 
source at location (0, 0, H) for which the source 
term takes the form 

Six,t)^Q6ix)S{y)S{z~H), (2) 

where Q is a constant emission rate [kgs^^] and S 
represents the Dirac delta function [m~^]. The co- 
ordinates X — (x, y, z) are chosen so that the x-axis 
is aligned parallel to the wind velocity. Although the 



wind velocity lies within the horizontal plane by our 
earlier assumption there is nonetheless a small 
vertical velocity component M^sct that corresponds 
to gravitational settling. As a result, the velocity ap- 
pearing in (III) takes the form u — {U{t), 0, — Wget), 
where U{t) is the wind speed and the constant set- 
tling velocity is determined for spherical particles 
using Stokes' law 

Wset = p.gdV(18M), (3) 

where 

p = particle density [kgm~'^], 

d — particle diameter [m] , 

= air viscosity — 1.8 x 10"^ [kgm~^ 

g = gravitational acceleration = 9.8 [ms^'^]. 
Diffusive transport in the x-direction is typically 
much smaller than convective transport by the wind; 
hence, the diffusion term involving may be ne- 
glected, and Ky (x) and (x) may be taken as func- 
tions of downwind distance only. As a result, Eq. ([T]) 
reduces to 
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dc , dc 



dc 

Wset^ 



+ Six,t). (4) 



dz 

^(k —\ + —(k — 

dy \ ^ dy J dz \ ^ dz ^ 

This equation could be modified to include vertical 
variations in both wind velocity and dispersion co- 
efficients following [27l28j , although such effects will 
not be considered in this paper. 

For the moment, we assume the domain is of infi- 
nite extent in the x and y directions, and the posi- 
tive z direction; consequently, we can apply the fol- 
lowing Dirichlet boundary conditions at infinity: 

c(x, ±00, z) — 0, (5) 
c(±oo, y, z) = 0, (6) 
c(x,y,oo) = 0. (7) 

The ground surface {z = 0) is where deposition oc- 
curs, and so we impose the following mixed (Robin 
type) condition on particle flux 
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where Wdcp > is the dry deposition velocity which 
is usually assumed constant and is determined from 
experiments. 

3.1. Ermak's solution 

The Eqs. (g])-® have been well-studied and Er- 
mak ^ derived the following analytical solution 
which is valid for x > 0: 
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X erfc I + \ , (9) 

where Wo = Wdcp — \Wset- The quantities ay^z{x) 
(having units of m) represent standard deviations of 
concentration in the y and z directions and can be 
expressed in terms of the eddy diffusivities as 
2 /"^ 

o'L(^) ^ ^ KyAx')dx'. 

This integral simplifies to i^T = a'^U /2x when the 
eddy diffusivities are constant, although in practice 



the a values need to be taken as functions of the 
downwind distance x in order to match observed 
plume concentrations. 

Suppose that the height of a receptor above the 
ground is z = h] then the key quantity of interest 
is the contaminant deposition flux at z = h which 
is given by the expression Wdcpc\z=h from Eq. ([5]), 
with the height replaced by z = h. It is essential to 
observe that the deposition flux depends linearly on 
the emission rate Q through Eq. ([9]) , a fact that will 
be exploited to great advantage in this work. 

3.2. Multiple sources and receptors 

Instead of just a single point source, suppose in- 
stead that we have a set of Ng sources, each with 
emission rate Qs [kgs^^] and location for s = 
1,2, Ns- Likewise, there are Nr receptors located 
at positions rj^., each with total accumulated depo- 
sition Dr [kg], for r = 1,2,..., Nr. On the Trail 
site depicted in Figure [U there are Ng = A sources 
and Nr = 9 receptors that are of interest to this 
study. In the derivation that follows, we note that 
there are many similarities with Calder's [7] gen- 
eral model framework for multiple-source emissions, 
except that he was unconcerned with deposition or 
any specific Gaussian plume approximation for the 
advection-diffusion equation. 

In order to apply the solution from ^ for a uni- 
directional wind, a new set of transformed coordi- 
nates must be defined for each source s that trans- 
late the source location to the origin and then rotate 
coordinates so that the transformed a;-axis is aligned 
parallel with the wind velocity. To this end, we de- 
fine new coordinates x' which are related to x via 
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(10) 



where 6 corresponds to the angle the wind direction 
vector makes with the x-axis (measured counter- 
clockwise) and TZ^e represents the matrix that ro- 
tates vectors through an angle —9 in the x, y-plane 
(see Figure[2|) . The resulting coordinates have an X-- 
axis that is aligned with the wind direction. It is con- 
venient to rewrite the concentration from ([9|) by fac- 
toring out the source emission rate and making the 
dependence on the source location and wind explicit, 
letting c{Xg) = Qsp{x;^s,d,U). Consequently, the 
deposition flux may be written as 

Wdcpc\z=h = WdcpQsP{x;^g,9,U)\z=h, 

and we may then calculate the total accumulation of 
contaminant in kg within a dustfall jar at location 
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^ X 

Fig. 2. Relationship between the original coordinates x and 
transformed coordinates x' for a given wind field having 
magnitude U and direction angle 9. 

X over a time interval of length At as 

D{xi) = ^(t^dcp QsAAt)p{x; ^,,9, U), 

s=l 

where A is the cross-sectional area A of the jar open- 
ing. With a slight change in notation, we can write 
the total mass of a single contaminant deposited at 
receptor location r]^ as 

Dr = WdopAAtY,Qsp{Vr; ^s,0, C/)U=ft„ (11) 
s=l 

which clearly indicates that the deposition is a linear 
combination of the emission rates Qs- 

3.3. Time-varying wind 

Eqs. (fTT|) are derived for a wind that is constant 
over the time interval of length At. In practice, the 
wind velocity varies with time and these variations 
have a significant impact on the plume dispersion 
and consequently also the deposition. The wind 
speed and direction, U{t) and 0{t), are therefore 
both functions of time and so the coordinate trans- 
form (jlOp is also time-dependent. In order to make 
use of the Ermak solution, we divide time into Nt 
equally-spaced intervals of length At = T/Nt, and 
assume that U{t) and 9{t) can both be approx- 
imated by piecewise constant functions on each 
interval. The effects of time variation may then 
be incorporated by solving a sequence of Gaussian 
plume problems of the form ^ and the total mass 
deposited at receptor r is calculated by summing 
the results over each time interval. 
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(12) 

where p{x; 6*", J7") represents the expression for 
the concentration, using averaged values of [/" and 
0" corresponding to the rt*'* time interval. 

The Gaussian plume model assumes steady-state 
conditions, and so it is important to emphasize that 
we are making a major assumption when applying 
it to problems in which the wind velocity varies with 
time. As a result, the time step At is a critical pa- 
rameter: it must be chosen large enough that the 
plume can be considered sufficiently developed that 
it has reached a steady state, and yet not too large 
that the wind is seriously under-resolved. We choose 
a value of At = 10 min, motivated by a number of 
considerations: 

- Since "typical" wind speeds lie between 1 and 
5 ms~^, a given plume will be advected over a 
distance of 600 to 3000 m, which is close to the 
1000 m maximum separation between any indi- 
vidual source-receptor pair. 

- Beychok |23^ advocates an averaging interval of 
10 minutes or less for the Gaussian plume model 
based on the observation that one-hour intervals 
can lead to over-prediction of the concentration 
by as much as a factor 2.5. 

- Hanna et al. [26 note that Ermak's solution is 
derived based on the assumption that all variables 
are averaged over time periods of approximately 
10 min in length. 

Wind data is available down to a resolution of min- 
utes, and so presents no limitation on the choice of 
At. 



3.4. Parameter values and wind data 

The physical parameters corresponding to the 
contaminant particles are summarized in Table [TJ 
Because both zinc and strontium are deposited 
as sulphates, the parameters actually correspond 
to ZnS04 and SrS04. While zinc is the primary 
element of interest in this study, we will see in Sec- 
tion|4]that strontium is an important tracer element 
that plays a useful role in determining the solution 
to the inverse problem. The values for particle di- 
ameter and deposition velocity are consistent with 
data from [29 30 31j, while the settling velocity is 
computed from Stokes' law ([3|). 
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Table 1 

Values of the physical parameters for the two contaminants of interest, where zinc and strontium actually appear in the form 
of sulphates - ZnS04 and SrS04. Note that the particle diameter and deposition velocity for Sr (marked *) have been set 
equal to those for Zn, in the absence of other more reliable estimates. 



Parameter 


Symbol 


Units 


Value for 

Zn 


species q = 

Sr 


Density 




kg m~^ 


3540 


3960 


Molar mass 




kg mol~^ 


0.161 


0.184 


Diameter 






5 


5* 


Deposition velocity 


dcp 


ms~^ 


0.005 


0.005* 


Settling velocity 
(from Eq. lO) 




ms^^ 


0.0027 


0.0030 



The heights of the four contaminant sources, 
corrected for plume rise, are Hs = [15, 35, 15, 15] 
while the nine receptors are located at heights 
hr = [0,10, 10,1,15,2,3,12,12]. The dustfaU jars 
are glass containers in the shape of circular cylin- 
ders having a diameter of 0.162 m, and so the area 
parameter used in the deposition calculation in 
Eq. (HID is .4 = 0.0206 m^. 

The other key parameters appearing in the model 
are the standard deviations ay^z which we assume 
take the form a{x) = ax{l + hx)^'^ customarily at- 
tributed to Briggs [32^ • The constants appearing in 
this expression depend on the atmospheric stability 
class, and if we assume class D as discussed earlier 
then the values of the parameters are j33j : 



for (Ty-. a = 0.08, b = 0.0001, c = 0.5, 
forcr^: a = 0.06, 6 = 0.0015, c = 0.5. 



The wind data for all simulations come from ac- 
tual meteorological measurements taken on site and 
are specified as point measurements at 10-minute in- 
tervals. Figure [3] depicts the typical distributions of 
wind direction and velocity for a representative one- 
month period. The bimodal nature of the wind dis- 
tribution in a direction aligned roughly parallel to 
the Columbia River valley is evident from the wind 
rose diagram. A significant portion of wind measure- 
ments exhibit a zero velocity which cannot be used 
directly in Eq. ([9]). It is therefore usual to introduce 
a "cut-off" velocity C/min such that whenever the 
wind satisfies U ^ f/min, no contribution is made to 
the deposition during that time interval. We employ 
a cut-off of [/min = 0.1 which limits the number of 
neglected time intervals to 10% of the total. 



N 



Max. wind=7.47 ni/s 




Wind speed histogram 




Speed (iTi/s) 



Fig. 3. Measured wind data over the one-month period June 
3-July 2, 2002. Top: Wind rose diagram, with the proportion 
of calm winds identified in the central circle. Bottom: Wind 
speed histogram. 

3.5. Sample forward computation 

Using the parameter values given in the previous 
section, we next calculate a typical distribution of 
the zinc ground-level concentrations using the time- 
varying plume solution (|12p . The emission rates for 
the four sources were taken to be Q = [35, 80, 5, 5] x 
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10'^ kgyr~^, and measured wind data at 10 minute 
intervals was used for the month of June 2002. The 
mass of zinc deposited at each receptor was then 
computed to 3 significant digits as 

£) =[ 17.4, 31.0, 15.7, 1.63, 21.7, 

(13) 

6.33, 2.54, 5.02, 7.74 ] x 10"^ kg. ^ ^ 

We also computed ground-level zinc concentrations 
on a 100 X 100 grid of points covering the entire 
Trail site, and the results are displayed as a contour 
plot in Figure ID The concentration clearly peaks 
at locations that lie close to sources SI and S2 and 
aligned with the prevailing wind direction. This is to 
be expected since SI and S2 have by far the largest 
emission rates among the four sources. 



Contours of Zn concentration (g/m ) 




500 1000 

X (m) 



Fig. 4. Forward simulation of the monthly cumulative de- 
position of zinc (in gm"'^) at each point in the domain for 
source emission rates Q = [35,80,5,5] X 10'^ kgyr~^. The 
wind data used to generate this plot corresponds to June 
2002, for which the dominant wind direction is towards the 
southeast. Note that compass north is directed vertically in 
this plot, which is rotated relative to the aerial photo in 
Fig. [I 

4. Inverse problem 

We now consider the inverse problem for which 
values of the zinc deposition Dr are specified at each 
receptor r = 1, 2, . . . , Nr, and the unknowns are the 
source emission rates Qs for s — 1,2, . . . , Ng. It is 
convenient to rewrite Eqs. (|12p in the more compact 
form 

D = FQ (14) 



where D and Q are vectors containing the deposi- 
tions and emission rates respectively, and P is an 
Nr X Ns matrix whose r, s entry is given by 

Nt 

Prs = WdcpAAt ^p(,7,;C„r,C/")U^;,^. (15) 

n=l 

Since there will typically be more receptor measure- 
ments than sources {Nr > Ng), this is clearly an 
overdetermined system of equations for Q; hence, 
there is not a unique solution and we can hope at 
best to obtain an approximation of Q. 

We employ a linear least squares method to de- 
termine a solution of Eq. (|14p . and in all of our 
computations we employ the Isqlin solver in MAT- 
LAB. A typical simulation requires only about 20 s 
of CPU time on a MacBook Pro with a 2.5 GHz Intel 
dual-core processor, and so our method is very effi- 
cient and well-suited to performing parametric stud- 
ies. MacKay et al. P3] used a similar least squares 
method to estimate parameters such as surface mass 
transfer rate and Peclet number rather than emis- 
sion rates. 

As a consistency check on our inverse algorithm, 
we consider emissions of a single contaminant (zinc) 
under the influence of wind data during the month of 
June 2002. Deposition values are taken equal to the 
outputs from the forward computation in Eq. (jl3p , 
with all other parameters chosen the same. As ex- 
pected, the emission rates calculated with the in- 
verse algorithm are identical to within round-off er- 
ror. 

4.1. Multiple contaminants and linear constraints 

As mentioned earlier, dustfall measurements are 
made of several contaminants, not just zinc. The 
contents of each receptor underwent a trace analy- 
sis that measures the mass of a number of trace ele- 
ments, of which we are interested in zinc, strontium 
and sulphur (in the form of SO4). Although zinc is 
the contaminant of primary interest in this study, 
the sulphur and strontium arise from the same chem- 
ical processes that produce zinc, and consequently 
we have included them in the inversion procedure 
with the hope that they will increase the accuracy 
of the zinc estimate. Because zinc and strontium are 
deposited as sulphates, the emission rate of sulphur 
cannot be included as an independent variable and 
so we need to estimate emission rates for only zinc 
and strontium. 
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By modifying our previous notation slightly, we 
can identify contaminant species with a superscript 
q in the emission rates Q1 (where q = Zn or Sr) and 
depositions -D^ {q = Zn, Sr or SO4). Our task is then 
to solve the following three linear systems 



Jun 2002 



D 



Sr 



D 



SO4 



pZn gZr 
pSr gSr 



,Zn 



SO4 



(16) 



Sr/-,Sr 



MS' 



where the matrices P' are given in psp and depend 
on the contaminant q through the deposition veloc- 
ity W^dep and settling velocity = g p" (^9)7X8^. 
Since the scenario we are studying has Ng = 4 
sources and Nr = 9 receptors, Eqs. (fT6|) represent 
an overdetermined system of 27 equations in 8 un- 
knowns. Once Q^" and Q^'^ are known, the emission 
rate for sulphate can be obtained from 



Q 



SO4 



M 



SO4 



M 



SO4 



-Q 



Sr 



(17) 



In addition to the atmospheric transport pro- 
cesses embodied by the above equations, the chemi- 
cal processes generating the contaminants introduce 
a number of further constraints on the emission 
rates. To be more specific, we note that source 
SI corresponds to a collection of zinc purification 
stacks, S2 is a sulphide leach plant cooling tower, 
and S3 and S4 are identical electrolytic cooling 
stacks. As a result, we can make the following as- 
sumptions: 

- Sources 3 and 4 are identical: 

Ql-Ql = 0, for q = Zn, Sr. (18) 

- The mass ratio of zinc to strontium in sources 3 
and 4 is approximately 6000 to 1: 

Qf' - 6000Qf = for s = 3, 4. (19) 

- No strontium is emitted from sources 1 and 2: 

gSr ^ ^Sr 



= 0. 



(20) 



All emission rates must be non-negative: 



iSr 



^0 for s = 1,2,3,4. 



(21) 



When taken together, Eqs. (|18p - ((2T|) represent 6 
equality and 8 inequality constraints that supple- 
ment the overdetermined system (jl6p . We can there- 
fore continue to make use of the MATLAB linear 
least squares solver Isqlin which conveniently per- 
mits the inclusion of both equality and inequality 
constraints. 
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Fig. 5. Experimental measurements of the mass of contam- 
inants accumulated within each receptor during the month 
of June 2002. The zinc and sulphate figures are shown in 
units of mg, while strontium is scaled by 10^ so that it is 
visible on the same axes. 

5. Numerical simulations 



5.1. Base case 

We begin by focusing on a number of inverse cal- 
culations using a single month of wind and depo- 
sition data corresponding to June 2002 - we refer 
to this simulation as the "base case." The receptor 
measurements used as inputs are displayed in Fig.[5l 
while the wind data is the same as that shown ear- 
lier in Fig. [H The estimated source emission rates 
are displayed in Fig. [6] in units of tonnes per year. 
We are primarily interested in zinc emissions, which 
come to a total of 92.44 Tyr~^ for all four sources. 
This total should be compared with the figure of 
116.48 T yr^^ reported by Teck-Cominco in the Na- 
tional Pollutant Release Inventory for the year 2002 
[25] . The fact that these two figures are so close gives 
us some confidence that our inverse approach for 
estimating emissions is a reasonable one. We have 
not attempted a comparison between the other two 
contaminants because strontium emissions are not 
reported to NPRI, and our model does not account 
for all sources of atmospheric sulphate emissions on 
the Trail site. 

5.2. Parameter sensitivities 

In this section, we investigate the sensitivity of 
our inverse solution to changes in several key pa- 
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Fig. 6. Estimated emissions rates for the "base case" cor- 
responding to June 2002 (the strontium value is scaled by 
10*). 

rameters, namely receptor height (/ir), source height 
(Hs) and deposition velocity (W^^^). These param- 
eters are singled out because they are all subject to 
significant variations for the following reasons: 

- There is a great deal of freedom in the placement 
of individual receptors, which for example can be 
located at ground level or else on the top of build- 
ings or other structure. 

- The height of each source must be adjusted to take 
into account the effect of plume rise, which can be 
estimated using empirical formulas but there re- 
mains a significant degree of uncertainty in these 
plume rise estimates. 

- Particle deposition velocities are known to de- 
pend strongly on atmospheric conditions and on 
whether deposition happens under dry or wet con- 
ditions. 

We perform a series of simulations, modifying each 
parameter by ±10% and ±20% from the base case 
value. The resulting emission estimates are com- 
pared in Figs. [21 [8] and [9] for parameters hr, Hs and 
^dcp respectively. Only zinc emission rates are de- 
picted here since similar levels of sensitivity are ex- 
perienced for the other contaminants. 

From these results, we conclude that the solu- 
tion is most sensitive to source height, where rela- 
tive changes on the order of 40-50% are seen. Much 
less variation is seen in response to changes in H 
and Wdcp- These sensitivities are consistent with 
Miller and Hively [31] , who reviewed a wide range of 
Gaussian plume type models and found that for ele- 
vated sources, ground-level concentrations are typi- 
cally estimated to within a factor of 0.65 to 1.35 of 



Fig. 7. Effect of cfianges in receptor height hr on Zn emission 
rates. 
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Fig. 8. Effect of changes in source height Hs on Zn emission 
rates. 
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Fig. 9. Effect of changes in deposition velocities VFj^p on Zn 
emission rates. 
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Fig. 10. Effect of random noise in deposition measurements 
on Zn emission rates. 



the actual values. 



5.4. A note on ill- conditioning 

Our approach does not suffer from the extreme 
levels of sensitivity observed by others in inverse 
computations based on Gaussian plume type solu- 
tions (for example, [34], [35] and [36]). We can ex- 
plain this apparent discrepancy by noting that our 
simulations are performed over relatively short time 
periods as well as being restricted to ground level 
and to areas very close to the emission source. In 
contrast, Enting and Newsam's study of the inverse 
emissions problem indicated that computed concen- 
trations at distances further than 10 km downwind 
from a source are highly sensitive when high-altitude 
pollutant measurements are used. Although these 
and other studies of atmospheric dispersion focus on 
long-range transport (over lO's or even lOO's of km) 
our approach benefits from the fact that we consider 
transport over much shorter spatial scales. 



5.3. Noise in deposition measurements 

Even though we expect that the emission rates 
should remain approximately constant over time (at 
least during a given year) there is nonetheless a 
significant spread in deposition measurements from 
month to month. There are a number of possible ex- 
planations for this variation, including measurement 
errors, wet deposition during rainy periods (where 
"scrubbing" can significantly reduce the amount de- 
posited) or contamination from secondary sources 
not accounted for in the model. Indeed, Goodarzi 
et al. [21] performed an experimental study of emis- 
sions at the Teck Cominco's Trail operation and de- 
termined that secondary sources on the site (such 
as ore concentrate and slag storage piles) experience 
resuspension of particles which may interfere with 
deposition measurements. 

In any case, it is important to understand the ef- 
fect that possible errors in particulate measurements 
may have on emission estimates. To this end, we 
take each deposition measurement and scale it by a 
normally distributed random number chosen from 
the interval [1 — a, 1 -I- a] for values of a = 0.1, 0.2 
and 0.3. The results are summarized in Fig. llOi from 
which it is clear that even for the largest noise ratio 
a = 0.3, the corresponding errors in the computed 
emission rates are relatively small. One possible ex- 
planation for this limited sensitivity to noise is that 
positive and negative contributions to errors in the 
input data may cancel each other out on average. 



5.5. Deposition over several months 

Deposition measurements are available for a total 
of six monthly periods: June, October and Novem- 
ber in 2001; and May, June and July in 2002. All 
of these measurements correspond to periods dur- 
ing which the Trail smelter was in continuous opera- 
tion without being shut down; therefore, one would 
expect the depositions to be fairly consistent from 
month to month. The measured values for zinc are 
summarized in Fig. 111! from which we observe that 
there are significant variations at each receptor loca- 
tion over time, although the trend is fairly consistent 
between one receptor and another. Furthermore, the 
largest accumulations are consistently measured in 
receptors located closest to the primary sources SI 
and S2, as expected. Similar trends are observed in 
the strontium and sulphate data and so they are not 
depicted. 

Keeping in mind our earlier discussion in Sec- 
tion 15.31 of the solution sensitivity to the error in 
deposition data, we next apply our inverse algo- 
rithm for each one-month periods individually. The 
resulting monthly emission rate estimates for zinc 
are given in the bar plot in Fig. 1121 There is con- 
siderable variation between the results for individ- 
ual months, which is not surprising considering the 
variations in the input data. Furthermore, two of the 
S2 estimates (corresponding to June 2001 and July 
2002) are identically zero, indicating that the least 
squares algorithm has forced the S2 value up against 
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Fig. 11. Simulated Zn deposition rates for every monthly pe- 
riod under consideration. For each receptor, the axis limits 
are drawn as horizontal dotted lines and are scaled so that 
the lower limit is zero and the upper limit represents the 
maximum deposition rate (given by the number in paren- 
theses to the left of each curve in Tyr~^). 
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Fig. 12. Computed Zn emission rates for each monthly pe- 
riod. 

the boundary of the inequality constraint Qf" ^ 0. 
These two results are clearly questionable because 
S2 is by far the largest source of zinc on the site. 

With an aim to explaining this discrepancy, we 
take a closer look at the zinc deposition data in 
Fig. [TT] and observe that the R3 deposition mea- 
surements are unusually high in relation to other 
nearby receptors (in particular, see the peak value of 
135 mg measured in June 2001). Instead, one would 
expect that the values at R3 and R4 are much closer 
to each other since their location relative to the pri- 
mary sources SI and S2 is so similar. Furthermore, 
it seems reasonable that the maximum deposition 
should occur at a point lying to the southeast or 
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Fig. 13. Computed Zn emission rates for each monthly pe- 
riod, leaving out receptor R3 (compare to Fig. I12| l. 



northwest of sources SI and S2, rather than at R3 
which is perpendicular to the direction of the pre- 
vailing winds. We conclude therefore that the large 
deposition at R3 is most likely attributable to mea- 
surement errors or some other anomaly, and hence 
R3 should be excluded from the calculations. 

Upon repeating the previous set of simulations 
with R3 excluded, we obtain the results shown in 
Fig. [131 All SI and S2 emission estimates are now 
nonzero and the anomalously high SI estimate for 
June 2002 is reduced to be in more line with the 
other values. Although there remain significant 
variations between a few of the monthly estimates 
(namely November 2001 and July 2002), these re- 
sults are much more reasonable. Since the primary 
quantity of interest in this study is the total emis- 
sions from all sources, another way of viewing these 
results is via the total emissions from all four sources 
for each monthly period, as shown in Fig. 1141 

If we assume that emissions are constant through- 
out the year, then we may also use all months of de- 
position data from a given year to estimate a single 
"aggregate" emission rate. The aggregate estimates 
based on the 2001 and 2002 data are shown in Ta- 
ble [2] along with the total emission rates for the in- 
dividual months. One advantage of using the aggre- 
gate calculation is that it serves to average out some 
of the variation between the individual monthly re- 
sults. The aggregate figures can then be compared to 
the values from the NPRI database of 100.33 T yr"! 
in 2001 and 116.48 Tyr^^ in 2002. For both years, 
the aggregate estimates are within 10-20% of the re- 
ported values, and an increasing trend from 2001 to 
2002 is also captured. 
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Fig. 14. Total emissions for each contaminant. The solid lines 
correspond to simulations with R3 included in the inverse 
calculation, while the dashed lines are without R3. 

Table 2 



Summary of total zinc emissions for each month of interest, 
the corresponding aggregate totals for each year, and the 
emissions reported to NPRI. Receptor R3 is omitted. 





Zn Emission rates (Tyr ^) 




2001 


2002 




130.5 


121.7 


Monthly estimates 


104.9 


112.1 




63.64 


40.97 


Aggregate estimates 


95.77 


104.5 


Reported to NPRlEl] 


100.33 


116.48 



6. Conclusions 

We have derived a linear least squares approach 
for estimating contaminant emissions from several 
point sources given time-varying wind data and 
monthly-averaged deposition measurements. The 
model is based on Ermak's Gaussian plume type 
solution to the advection-diffusion equations [S] 
which incorporates particle settling and surface de- 
position. The novelty of our approach stems from 
its focus on short range emissions (within 1000 m of 
the source), incorporating additional equality and 
inequality constraints on contaminant sources, and 
its ability to handle time- varying wind data. 

The algorithm is used to estimate emissions of 
several contaminants from a smelting operation in 
Trail, British Columbia using dustfall measurements 
taken over several one-month periods during 2001- 
2002. We demonstrate that the algorithm is efhcient 
and robust to changes in several key parameter val- 
ues in comparison to other approaches published in 



the literature. While there is significant variation 
from month-to-month and between the individual 
sources, the average annual estimate of the emission 
rate is within 10-20% of the values reported to En- 
vironment Canada for the years in question. 

There remain a number of aspects of the model 
which require further study and refinement, for ex- 
ample by using more careful estimates of the contri- 
bution of plume rise to source heights, incorporating 
the effects of wet deposition during rainy periods, 
and allowing further variation in atmospheric con- 
ditions through the Pasquill-Gifford stability classi- 
fication. We would also like to apply our insight into 
contaminant plume dynamics to guide the design of 
future studies of the Trail site; in particular, by re- 
locating a receptor such as R3 to a more advanta- 
geous position (e.g., within the areas of peak con- 
centration southeast or northwest of sources SI and 
S2). Finally, we will perform a direct finite volume 
discretization of the 3D advection-diffusion equa- 
tion, for which reliable estimates for the turbulent 
eddy diffusivities are essential. Comparisons of these 
simulations with the results of our inverse Gaussian 
plume algorithm will appear in a companion publi- 
cation [37] . 
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